{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "988c4ea4",
   "metadata": {},
   "source": [
    "# Film preperation using particle deposition.\n",
    "\n",
    "This tutorial demonstrates particle film generation using ballistic deposition.\n",
    "\n",
    "#### Current maintainer\n",
    "- Stefan Endres (s.endres@iwt.uni-bremen.de)\n",
    "\n",
    "#### Previous authors:     \n",
    "- Valentine Baric \n",
    "- Lutz Mädler$^*$ (lmaeder@iwt.uni-bremen.de)\n",
    "- Norbert Riefler \n",
    "\n",
    "#### Literature: \n",
    "##### Film generation primary sources:\n",
    "\n",
    "- Mädler et al. (2006), Nanotechnology 17, 3783-4795\n",
    "- Norbert, R. and Mädler. L. (2010) Journal of Nanoparticle Research 12.3 : 853-863.\n",
    "\n",
    "##### Secondary:\n",
    "\n",
    "- Filippov et al. (2000), Journal of Colloid and Interface Science 229, 261-273\n",
    "- Ermak and Buckholz (1980), Journal of Computational Physics 35. 169-182\n",
    "- Chan and Dahneke (1981), J. Appl. Phys. 52 3106-10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "89834caa",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Imports\n",
    "# Std. library\n",
    "import sys, random, os, datetime, warnings, socket, argparse\n",
    "import multiprocessing as mp\n",
    "import numpy as np\n",
    "# Dependencies:\n",
    "from scipy import spatial, stats\n",
    "import shutil, fcntl, termios, struct\n",
    "# Local imports:\n",
    "from main_function import *\n",
    "import functions as f\n",
    "from create_aggregate import create_aggregate as SA\n",
    "from create_aggregate_CCA import create_aggregate as CCA\n",
    "import debugging as debug\n",
    "from classes_agg import Aggregate as Aggregate\n",
    "from classes_agg import Particle as Particle\n",
    "import write_output\n",
    "import terminal_print\n",
    "import statistics\n",
    "warnings.filterwarnings('error')\n",
    "if sys.stdout.isatty():\n",
    "    Redirection = 0\n",
    "    COLS = struct.unpack('hh',  fcntl.ioctl(sys.stdout, termios.TIOCGWINSZ, '1234'))[1] -1\n",
    "else:\n",
    "    COLS = 60\n",
    "    Redirection = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f18f6920",
   "metadata": {},
   "source": [
    "# Parameters\n",
    "\n",
    "Next we import the parameters for the film generation, the `dda_template.py` file can be modified to add custom parameters, alteratively the desired imports can be modified in this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b5c1c70f",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Exception ignored in: <function Pool.__del__ at 0x7f1e482603a0>\n",
      "Traceback (most recent call last):\n",
      "  File \"/usr/lib/python3.10/multiprocessing/pool.py\", line 265, in __del__\n",
      "    _warn(f\"unclosed running multiprocessing pool {self!r}\",\n",
      "ResourceWarning: unclosed running multiprocessing pool <multiprocessing.pool.Pool state=RUN pool_size=1>\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Aggregate 0/1 complete\n"
     ]
    }
   ],
   "source": [
    "# The inputfile has to be copied to the folder of the source, otherwise it cannot be imported (shutil.copy2) (only if they are not already in the same folder)\n",
    "import dda_template as params\n",
    "source_file = os.getcwd()\n",
    "PID = os.getpid() \n",
    "params_file = 'dda_template.py'\n",
    "\n",
    "# Initiate parameters:\n",
    "init_params(params_file, source_file)\n",
    "# Example parameters that can be modified:\n",
    "params.N_tot = 20#0000  # Maximum number of primary particles (PP)  (test with low number)\n",
    "params.experiment  # The name of this experiment\n",
    "params.output_folder  # Folder for output data\n",
    "params.T  # Gas temperature\n",
    "params.D_f  # Fractal dimension\n",
    "k_f = params.k_f  # Prefractor for D_f\n",
    "epsilon = params.epsilon  # Minimum distance between two PP\n",
    "delta = params.delta  # Maximum distance between two PP\n",
    "deformation = params.deformation  # Deformation factor\n",
    "params.rho  # Density of the particle\n",
    "params.dt  # Timestep\n",
    "Height_above = params.Height_above  # Initialization hight above the highest particle\n",
    "\n",
    "# Parrallelization:\n",
    "num_cores = 1  # Specify the number of CPUs used for parallel computing\n",
    "jobserver = mp.Pool(processes = num_cores)\n",
    "if num_cores == -1:\n",
    "    num_cores = mp.cpu_count()\n",
    "elif num_cores <= 0 or num_cores > mp.cpu_count():\n",
    "    num_cores = mp.cpu_count()\n",
    "    print(\"\"\"Invalid number of CPUs\n",
    "        ==> Maximum number of CPU used ('%i')\"\"\" %num_cores)\n",
    "    \n",
    "jobserver = mp.Pool(processes = num_cores)\n",
    "\n",
    "# Algorithmic controls:\n",
    "k_B = 1.3806488e-23         # Boltzmann-Constant in J/K\n",
    "start = datetime.datetime.now().replace(microsecond=0) # Starting time\n",
    "cwd = os.getcwd()           # Current directory\n",
    "RESTART = False\n",
    "random.seed(params.custom_seed_traj)  # set a new seed for random numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "795d89d2",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Exception ignored in: <_io.FileIO name='./dda/log.dda' mode='wb' closefd=True>\n",
      "Traceback (most recent call last):\n",
      "  File \"/tmp/ipykernel_614345/2376010603.py\", line 4, in <cell line: 4>\n",
      "ResourceWarning: unclosed file <_io.TextIOWrapper name='./dda/log.dda' mode='w' encoding='UTF-8'>\n"
     ]
    }
   ],
   "source": [
    "archive_folder = './archive'\n",
    "\n",
    "#current_date = now.strftime(\"%Y%m%d\")\n",
    "write_lammpstrj = output_trj(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fbe0b20",
   "metadata": {},
   "source": [
    "# Film generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bdd20042",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "************************************************************\n",
      "Process data\n",
      "\tExperiment:\t\tdda\n",
      "\tStarting Time:\t\t2022-06-16 17:50\n",
      "\tPID:\t\t\t614345\n",
      "\tCPUs:\t\t\t1\n",
      "\n",
      "************************************************************\n",
      "System parameters\n",
      "\tPe:\t\t\t1.00e+00\n",
      "\tBox size:\t\t280\n",
      "\tPrimary Particles:\t20\n",
      "\tCluster per Aggregate:\t6\n",
      "\tAggregates:\t\t1\n",
      "\tParticle Number Dist.:\tnone\n",
      "\tParticle Size Dist.:\tlognorm\n",
      "\tMedian radius (nm):\t4.479e+00\n",
      "\tSigma:\t\t\t3.700e-01\n",
      "\tNeighbor Mode:\t\t2\n",
      "\tFriction Model:\tChan&Dahneke\n",
      "\n",
      "************************************************************\n",
      "Modules\n",
      "\tCalculate Coordination (2)\n",
      "\t\tThreshold: Delta (2)\n",
      "\tCalculate Percolation (3)\n",
      "\t\tThreshold: One rp (1)\n",
      "\n",
      "************************************************************\n",
      "Outputfiles\n",
      "\tDefault outputdata\n",
      "\tlog file (1)\n",
      "\tLIGGGHTS hybrid granular/molecular data file (2)\n",
      "\tLIGGGHTS granular data file (3)\n",
      "\tVTK Output for paraview (5)\n",
      "\tHistogram with primary particle size distibution (7)\n",
      "\tHistogram with primary particle number distibution (8)\n",
      "\n",
      "************************************************************\n",
      "Debugging routines:\n",
      "\tCalculate overlapping particles with boundaries (4)\n",
      "\n",
      "************************************************************\n",
      "\n",
      "##### Aggregate generation initiated\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "IOStream.flush timed out\n"
     ]
    }
   ],
   "source": [
    "results, aggs, N_PP_per_Agg = generate_film(params, source_file, num_cores, jobserver)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "72de6706",
   "metadata": {},
   "outputs": [],
   "source": [
    "pt, maximal_radius = global_pp_indexing(aggs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "19f02976",
   "metadata": {},
   "outputs": [],
   "source": [
    "write_output.N_histogram(N_PP_per_Agg, params.output_folder, params.experiment, params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a5755ab",
   "metadata": {},
   "source": [
    "# Particle Transport\n",
    "At first the velocity is calculated according to the Langevin\n",
    "equation of motion. This equation is solved according to Ermak and\n",
    "Buckholz, 1980 with a random fluctuating force (B1/B2) and an applied force\n",
    "(only relevant for the z-direction (force_z)). According to this velocity\n",
    "the displacement is calculated and the aggregate is moved. Each timestep (v)\n",
    "is computed on basis of the previous step (v0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d4c0e216",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "##### Particle Transport calculation initiated\n",
      "\n"
     ]
    }
   ],
   "source": [
    "max_z = particle_transport(aggs, pt, params, maximal_radius, write_lammpstrj)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62228952",
   "metadata": {},
   "source": [
    "# Post processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a1ed4fce",
   "metadata": {},
   "outputs": [],
   "source": [
    "write_lammpstrj.close()  # Note, if you want to re-run the cells above this needs to be re-opened\n",
    "\n",
    "# Generate .trj and .vtk files in the output_folder    \n",
    "post_proc(params, pt, aggs, source_file, max_z)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be491f48",
   "metadata": {},
   "source": [
    "# Film properties and visualisations\n",
    "\n",
    "In the `functions.py` file functions for computing the film `packing_density`, `coordination` number and `percolation` of the film are found. The `.trj` and `.vtk` files can be visualised. The cell below allows for visualisation of the particle film using polyscope which can be downloaded easily from PyPi using `$ pip install polyscope`:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ecf4b9ab",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "### Reading Particle Structure File\n",
      "# Filename: ./dda/dda.trj\n",
      "    File ./dda/dda.trj using units nano and type dda\n",
      "# File Format: trj\n",
      "# Units: nano\n",
      "# Primary Particles: 20\n",
      "# Aggregates: 0\n",
      "# Box Width: 1.254e+03 x 1.254e+03\n",
      "# Film Height: 8.936e+01\n",
      "### Done\n",
      "\n",
      "\n",
      "[polyscope] Backend: openGL3_glfw -- Loaded openGL version: 4.6 (Core Profile) Mesa 22.1.1\n"
     ]
    }
   ],
   "source": [
    "from classes.particlelayer import *\n",
    "import numpy as np\n",
    "import json\n",
    "import polyscope as ps\n",
    "\n",
    "file  = './dda/dda.trj'\n",
    "layer = particlelayer(file, nodescription=False).layer\n",
    "\n",
    "# Shape\n",
    "layer._data.shape\n",
    "\n",
    "# Extract positions in space\n",
    "points = np.zeros([layer.data.shape[0], 3])\n",
    "points[:, 0] = layer.data[:, layer.lib['x']]\n",
    "points[:, 1] = layer.data[:, layer.lib['y']]\n",
    "points[:, 2] = layer.data[:, layer.lib['z']]\n",
    "\n",
    "# Extract the radius\n",
    "radii  = np.zeros(layer.data.shape[0])\n",
    "radii = layer.data[:, layer.lib['r']]\n",
    "\n",
    "# Plot in polyscope\n",
    "ps.init()\n",
    "ps.set_up_dir(\"z_up\")\n",
    "ps_cloud = ps.register_point_cloud(\"Particles\", points)\n",
    "ps_cloud.add_scalar_quantity(\"Particle radii\", radii)\n",
    "ps_cloud.set_point_radius_quantity(\"Particle radii\", autoscale=False)\n",
    "ps_cloud.set_color((1.0, 1.0, 1.0))\n",
    "ps.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f231224",
   "metadata": {},
   "source": [
    "Explanation of main classes and functions:\n",
    "#### AGGREGATE-CLASS\n",
    "    \n",
    "    Additional to this main program there are some more module files:\n",
    "    Whenever a property of an aggregate is needed you will see something like:\n",
    "    aggs[k].r\n",
    "    This means that in aggs there are all aggregates stored and \"k\" is the\n",
    "    index of the aggregate. \".r\" is the property of the aggregate, in this case\n",
    "    r shows all the radii of the PP within this aggregate.\n",
    "    All properties and functions which start like this (aggs[k].) are found in classes.py\n",
    "\n",
    "#### PARTICLE-CLASS\n",
    "    \n",
    "    Same as above for AGGREGATE. The abbreviation is pt[k]. It is also found in classes.py\n",
    "\n",
    "#### FUNCTIONS\n",
    "    \n",
    "    Functions which are not specific for one aggregate are in the document functions.py\n",
    "    So if there is a function like \"f.rank_z()\" you'll find it there."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f31482d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
